home *** CD-ROM | disk | FTP | other *** search
/ Nebula 1 / Nebula One.iso / Graphics / Plotting / aa_Intel_Only / VectorField.app / dipoleElectricField.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-01-02  |  3.6 KB  |  120 lines
  1. #include <stdio.h>
  2. #include <math.h>
  3. #define X1 0.46
  4. #define X2 0.54
  5. #define Y1 0.47
  6. #define Y2 0.53
  7. #define CHARGE1  0.1
  8. #define CHARGE2 -0.1
  9. #define SCALE .1
  10.  
  11. main()
  12. {
  13.     double x,y,z,u,v;
  14.   double x1=X1,x2=X2, y1=Y1,y2=Y2, q1=CHARGE1,q2=CHARGE2;
  15.     double dist1,dist2,E1x,E1y,E2x,E2y;
  16.     double maguv,scale,ratio, rmin=1.0,rmax = 0.0;
  17.     char title[81];
  18.     char c;
  19.     FILE *datfile;
  20.     
  21.     datfile=fopen("dipoleElectricField.data","w");
  22.     
  23.     for (z = 0.0; z < 0.201; z+=0.02){
  24.             /* title line */
  25.         fprintf(datfile,"z = %f, q1 = %f, q2 = %f\n", z,q1,q2);
  26.             /* charge centers (with zero-length vectors) */
  27.         fprintf(datfile,"%12.7f%12.7f%12.7f%12.7f%c\n", x1,y1,0.,0.,'r'); //red
  28.         fprintf(datfile,"%12.7f%12.7f%12.7f%12.7f%c\n", x2,y2,0.,0.,'b'); //blue
  29.         for (x = 0.1; x < .95; x+=0.2){
  30.             for (y = 0.1; y < .95; y+=0.2){
  31.                 if (y > y1 && y < y2 && x > x1 && x < x2) {
  32.                     /* do nothing, to avoid large vectors between the charges */
  33.                 } else {
  34.                     dist1=sqrt((x-x1)*(x-x1) + (y-y1)*(y-y1) + z*z);
  35.                     dist2=sqrt((x-x2)*(x-x2) + (y-y2)*(y-y2) + z*z);
  36.                     E1x = q1*(x-x1)/(dist1*dist1*dist1);
  37.                     E1y = q1*(y-y1)/(dist1*dist1*dist1);
  38.                     E2x = q2*(x-x2)/(dist2*dist2*dist2);
  39.                     E2y = q2*(y-y2)/(dist2*dist2*dist2);
  40.                     u=E1x+E2x;
  41.                     v=E1y+E2y;
  42.                     fprintf(datfile,"%12.7f%12.7f%12.7f%12.7f%c\n", x,y,u,v,'g');//green
  43.                 }
  44.             }
  45.         }
  46.         fprintf(datfile,"%12.7f%12.7f%12.7f%12.7f%c\n", 0.,0.,0.,0.,'q');   
  47.         //separator between plots for different z's
  48.     }
  49.     fclose(datfile);
  50.     
  51.     /* Now read values back to find rmin and rmax to scale vectors */
  52.  
  53.     datfile=fopen("dipoleElectricField.data","r");
  54.     
  55.     for (z = 0.0; z < 0.201; z+=0.02){
  56.         if (fgets(title, 81, datfile) == NULL) {
  57.             printf("Done (or error in reading title)\n");
  58.             break;
  59.         }
  60.         while (fscanf(datfile, "%lf%lf%lf%lf%c", &x,&y,&u,&v,&c)==5) {
  61.             /* note that x...v are doubles, hence %lf in fscanf */
  62.             //printf("x = %8.5f, y = %8.5f, u = %8.5f, v = %8.5f, c = %c\n",
  63.             //    x,y,u,v,c);
  64.             if (c == 'q') {
  65.                 break;
  66.             }
  67.             if (u == 0.0 && v == 0.0) {
  68.                 //do nothing, zero-length vector, will be ignored
  69.             } else {
  70.                 maguv = sqrt(u*u+v*v);
  71.                 if (maguv < rmin) rmin = maguv;
  72.                 if (maguv > rmax) rmax = maguv;
  73.             }
  74.         }
  75.         printf("Exit block with c = %c, rmin = %f, rmax = %f\n", c,rmin,rmax);
  76.         fscanf(datfile,"\n");
  77.     }     //end of for-loop on z
  78.     
  79.     printf("rmin = %f, rmax = %f\n", rmin,rmax);
  80.     scale = SCALE*log(rmax/rmin);
  81.     printf("scale = %f\n", scale);
  82.     fclose(datfile);
  83.     
  84.     /* Now re-calculate and write scaled values back into data file */
  85.  
  86.     datfile=fopen("dipoleElectricField.data","w");
  87.     for (z = 0.0; z < 0.201; z+=0.02){
  88.             /* title line */
  89.         fprintf(datfile,"z = %f, q1 = %f, q2 = %f\n", z,q1,q2);
  90.             /* charge centers (with zero-length vectors) */
  91.         fprintf(datfile,"%12.7f%12.7f%12.7f%12.7f%c\n", x1,y1,0.,0.,'r'); //red
  92.         fprintf(datfile,"%12.7f%12.7f%12.7f%12.7f%c\n", x2,y2,0.,0.,'b'); //blue
  93.         for (x = 0.1; x < .95; x+=0.2){
  94.             for (y = 0.1; y < .95; y+=0.2){
  95.                 if (y > y1 && y < y2 && x > x1 && x < x2) {
  96.                     /* do nothing, to avoid large vectors between the charges */
  97.                 } else {
  98.                     dist1=sqrt((x-x1)*(x-x1) + (y-y1)*(y-y1) + z*z);
  99.                     dist2=sqrt((x-x2)*(x-x2) + (y-y2)*(y-y2) + z*z);
  100.                     E1x = q1*(x-x1)/(dist1*dist1*dist1);
  101.                     E1y = q1*(y-y1)/(dist1*dist1*dist1);
  102.                     E2x = q2*(x-x2)/(dist2*dist2*dist2);
  103.                     E2y = q2*(y-y2)/(dist2*dist2*dist2);
  104.                     u = (E1x+E2x);
  105.                     v = (E1y+E2y);
  106.                     maguv = sqrt(u*u+v*v);
  107.                     ratio = SCALE*log(maguv/rmin)/log(rmax/rmin);
  108.                     u = u*ratio;
  109.                     v = v*ratio;
  110.  
  111.                     fprintf(datfile,"%12.7f%12.7f%12.7f%12.7f%c\n", x,y,u,v,'g');//green
  112.                 }
  113.             }
  114.         }
  115.         fprintf(datfile,"%12.7f%12.7f%12.7f%12.7f%c\n", 0.,0.,0.,0.,'q');   
  116.         //separator between plots for different z's
  117.     }
  118.     fclose(datfile);
  119. }
  120.